{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convex approximate dynamic programming\n",
    "\n",
    "We consider a stochastic control problem of the form\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & \\underset{T \\to \\infty}\\lim {\\mathbb E} \\left[\\frac{1}{T} \\sum_{t=0}^{T-1} \\|{x_t}\\|_2^2 + \\|{\\phi(x_t)}\\|_2^2\\right]\\\\[.2cm]\n",
    "\\mbox{subject to} & x_{t+1} = Ax_t + B\\phi(x_t) + \\omega_t,\n",
    "\\end{array}\n",
    "\\label{eq:adp}\n",
    "\\end{equation}\n",
    "where $x_t\\in\\mathbf{R}^n$ is the state, $\\phi:\\mathbf{R}^n \\to \\mathcal U \\subseteq \\mathbf{R}^m$ is\n",
    "the policy, $\\mathcal U$ is a convex set representing the allowed set of controls,\n",
    "and $\\omega_t\\in\\Omega$ is a (random, i.i.d.) disturbance.\n",
    "Here the variable is the policy $\\phi$, and the expectation is taken over\n",
    "disturbances and the initial state $x_0$. If $\\mathcal U$ is not an affine\n",
    "set, then this problem is in general very difficult to solve.\n",
    "\n",
    "A common heuristic for solving stochastic control problems is\n",
    "approximate dynamic programming (ADP), which parametrizes $\\phi$\n",
    "and replaces the minimization over functions $\\phi$ with a minimization over parameters.\n",
    "In this example, we take $\\mathcal U$ to be the unit ball and we represent $\\phi$\n",
    "as a particular quadratic *control-Lyapunov* policy.\n",
    "Evaluating $\\phi$ corresponds to solving the SOCP\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & u^T P u + x_t^T Q u + q^T u \\\\\n",
    "\\mbox{subject to} & \\|{u}\\|_2 \\leq 1,\n",
    "\\end{array}\n",
    "\\label{eq:policy}\n",
    "\\end{equation}\n",
    "with variable $u$ and parameters $P$, $Q$, $q$, and $x_t$. We can run\n",
    "gradient descent (SGD) on $P$, $Q$, and $q$ to\n",
    "approximately solve the original problem, which requires requires differentiating\n",
    "through the quadratic policy. Note that if $u$ were unconstrained, the original problem\n",
    "could be solved exactly, via linear quadratic regulator (LQR) theory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "import torch\n",
    "from scipy.linalg import solve_discrete_are\n",
    "from scipy.linalg import sqrtm\n",
    "\n",
    "from cvxpylayers.torch import CvxpyLayer\n",
    "\n",
    "# Generate data\n",
    "torch.manual_seed(1)\n",
    "np.random.seed(1)\n",
    "\n",
    "n = 2\n",
    "m = 3\n",
    "\n",
    "A = np.eye(n) + 1e-2 * np.random.randn(n, n)\n",
    "B = 1e-2 / 3 * np.random.randn(n, m)\n",
    "Q = np.eye(n)\n",
    "R = np.eye(m)\n",
    "\n",
    "# Compute LQR control policy\n",
    "P_lqr = solve_discrete_are(A, B, Q, R)\n",
    "P = R + B.T@P_lqr@B\n",
    "P_sqrt_lqr = sqrtm(P)\n",
    "\n",
    "# Construct CVXPY problem and layer\n",
    "x_cvxpy = cp.Parameter((n, 1))\n",
    "P_sqrt_cvxpy = cp.Parameter((m, m))\n",
    "P_21_cvxpy = cp.Parameter((n, m))\n",
    "q_cvxpy = cp.Parameter((m, 1))\n",
    "\n",
    "u_cvxpy = cp.Variable((m, 1))\n",
    "y_cvxpy = cp.Variable((n, 1))\n",
    "\n",
    "objective = .5 * cp.sum_squares(P_sqrt_cvxpy @ u_cvxpy) + x_cvxpy.T @ y_cvxpy + q_cvxpy.T @ u_cvxpy\n",
    "problem = cp.Problem(cp.Minimize(objective), [cp.norm(u_cvxpy) <= 1, y_cvxpy == P_21_cvxpy @ u_cvxpy])\n",
    "assert problem.is_dpp()\n",
    "policy = CvxpyLayer(\n",
    "    problem, [x_cvxpy, P_sqrt_cvxpy, P_21_cvxpy, q_cvxpy], [u_cvxpy])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we train the policy and plot the estimated average cost for each iteration of\n",
    "SGD for a numerical example, with $x \\in\n",
    "\\mathbf{R}^2$ and $u \\in \\mathbf{R}^3$, a time horizon of $T=25$, and a batch\n",
    "size of $8$. We initialize our policy's parameters with the LQR solution,\n",
    "ignoring the constraint on $u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(iter 0) loss: 1.51906 \n",
      "(iter 1) loss: 1.51596 \n",
      "(iter 2) loss: 1.51004 \n",
      "(iter 3) loss: 1.50076 \n",
      "(iter 4) loss: 1.48964 \n",
      "(iter 5) loss: 1.47732 \n",
      "(iter 6) loss: 1.46409 \n",
      "(iter 7) loss: 1.44985 \n",
      "(iter 8) loss: 1.43349 \n",
      "(iter 9) loss: 1.41393 \n",
      "(iter 10) loss: 1.39383 \n",
      "(iter 11) loss: 1.37273 \n",
      "(iter 12) loss: 1.35031 \n",
      "(iter 13) loss: 1.32644 \n",
      "(iter 14) loss: 1.2979 \n",
      "(iter 15) loss: 1.26749 \n",
      "(iter 16) loss: 1.23994 \n",
      "(iter 17) loss: 1.21162 \n",
      "(iter 18) loss: 1.18454 \n",
      "(iter 19) loss: 1.15832 \n",
      "(iter 20) loss: 1.13336 \n",
      "(iter 21) loss: 1.10918 \n",
      "(iter 22) loss: 1.08801 \n",
      "(iter 23) loss: 1.07002 \n",
      "(iter 24) loss: 1.05461 \n",
      "(iter 25) loss: 1.04056 \n",
      "(iter 26) loss: 1.02755 \n",
      "(iter 27) loss: 1.0159 \n",
      "(iter 28) loss: 1.00493 \n",
      "(iter 29) loss: 0.995569 \n",
      "(iter 30) loss: 0.987571 \n",
      "(iter 31) loss: 0.980465 \n",
      "(iter 32) loss: 0.974158 \n",
      "(iter 33) loss: 0.968241 \n",
      "(iter 34) loss: 0.962846 \n",
      "(iter 35) loss: 0.958006 \n",
      "(iter 36) loss: 0.953747 \n",
      "(iter 37) loss: 0.949987 \n",
      "(iter 38) loss: 0.946659 \n",
      "(iter 39) loss: 0.943703 \n",
      "(iter 40) loss: 0.94107 \n",
      "(iter 41) loss: 0.938717 \n",
      "(iter 42) loss: 0.936607 \n",
      "(iter 43) loss: 0.934709 \n",
      "(iter 44) loss: 0.932997 \n",
      "(iter 45) loss: 0.931448 \n",
      "(iter 46) loss: 0.930042 \n",
      "(iter 47) loss: 0.928762 \n",
      "(iter 48) loss: 0.927593 \n",
      "(iter 49) loss: 0.926524 \n",
      "(iter 50) loss: 0.925542 \n",
      "(iter 51) loss: 0.924639 \n",
      "(iter 52) loss: 0.923806 \n",
      "(iter 53) loss: 0.923036 \n",
      "(iter 54) loss: 0.922322 \n",
      "(iter 55) loss: 0.921659 \n",
      "(iter 56) loss: 0.921042 \n",
      "(iter 57) loss: 0.920465 \n",
      "(iter 58) loss: 0.919926 \n",
      "(iter 59) loss: 0.919421 \n",
      "(iter 60) loss: 0.918946 \n",
      "(iter 61) loss: 0.918499 \n",
      "(iter 62) loss: 0.918078 \n",
      "(iter 63) loss: 0.917679 \n",
      "(iter 64) loss: 0.917302 \n",
      "(iter 65) loss: 0.916944 \n",
      "(iter 66) loss: 0.916604 \n",
      "(iter 67) loss: 0.916281 \n",
      "(iter 68) loss: 0.915972 \n",
      "(iter 69) loss: 0.915678 \n",
      "(iter 70) loss: 0.915396 \n",
      "(iter 71) loss: 0.915126 \n",
      "(iter 72) loss: 0.914867 \n",
      "(iter 73) loss: 0.914618 \n",
      "(iter 74) loss: 0.91438 \n",
      "(iter 75) loss: 0.914149 \n",
      "(iter 76) loss: 0.913928 \n",
      "(iter 77) loss: 0.913714 \n",
      "(iter 78) loss: 0.913507 \n",
      "(iter 79) loss: 0.913307 \n",
      "(iter 80) loss: 0.913114 \n",
      "(iter 81) loss: 0.912927 \n",
      "(iter 82) loss: 0.912745 \n",
      "(iter 83) loss: 0.912569 \n",
      "(iter 84) loss: 0.912398 \n",
      "(iter 85) loss: 0.912232 \n",
      "(iter 86) loss: 0.91207 \n",
      "(iter 87) loss: 0.911913 \n",
      "(iter 88) loss: 0.91176 \n",
      "(iter 89) loss: 0.91161 \n",
      "(iter 90) loss: 0.911465 \n",
      "(iter 91) loss: 0.911323 \n",
      "(iter 92) loss: 0.911185 \n",
      "(iter 93) loss: 0.911049 \n",
      "(iter 94) loss: 0.910917 \n",
      "(iter 95) loss: 0.910788 \n",
      "(iter 96) loss: 0.910662 \n",
      "(iter 97) loss: 0.910539 \n",
      "(iter 98) loss: 0.910418 \n",
      "(iter 99) loss: 0.910299 \n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def train(iters):\n",
    "    # Initialize with LQR control lyapunov function\n",
    "    P_sqrt = torch.from_numpy(P_sqrt_lqr).requires_grad_(True)\n",
    "    P_21 = torch.from_numpy(A.T @ P_lqr @ B).requires_grad_(True)\n",
    "    q = torch.zeros((m, 1), dtype=torch.double, requires_grad=True)\n",
    "    variables = [P_sqrt, P_21, q]\n",
    "    A_tch, B_tch, Q_tch, R_tch = map(torch.from_numpy, [A, B, Q, R])\n",
    "\n",
    "    def g(x, u):\n",
    "        return (x.t() @ Q_tch @ x + u.t() @ R_tch @ u).squeeze()\n",
    "\n",
    "    def evaluate(x0, P_sqrt, P_21, q, T):\n",
    "        x = x0\n",
    "        cost = 0.\n",
    "        for _ in range(T):\n",
    "            u, = policy(x, P_sqrt, P_21, q)\n",
    "            cost += g(x, u) / T\n",
    "            x = A_tch @ x + B_tch @ u + .2 * torch.randn(n, 1).double()\n",
    "        return cost\n",
    "\n",
    "    def eval_loss(N=8, T=25):\n",
    "        return sum([evaluate(torch.zeros(n, 1).double(), P_sqrt, P_21, q, T=T)\n",
    "                    for _ in range(N)]) / N\n",
    "\n",
    "    results = []\n",
    "    optimizer = torch.optim.SGD(variables, lr=.02, momentum=.9)\n",
    "    for i in range(iters):\n",
    "        # use same seeds each iteration to get pretty training plot\n",
    "        torch.manual_seed(1)\n",
    "        np.random.seed(1)\n",
    "        optimizer.zero_grad()\n",
    "        loss = eval_loss()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        results.append(loss.item())\n",
    "        print(\"(iter %d) loss: %g \" % (i, results[-1]))\n",
    "    return results\n",
    "\n",
    "\n",
    "results = train(iters=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhcZZn38e/dVb131u7OnpA9IWFL6CABAkFAA8gOsrgxhokMOug7rzrqKOPyjsqgM44biixBBlBRWUUDghiWENJANpJAQsi+dGfrLJ3equ73jzoNndCdriR1+nR3/T7XVVfVWarOfTzYv5znnPM85u6IiEj2yom6ABERiZaCQEQkyykIRESynIJARCTLKQhERLJcPOoCDldZWZkPHz486jJERLqUV199dZu7l7e2rMsFwfDhw6msrIy6DBGRLsXM1ra1TE1DIiJZTkEgIpLlFAQiIllOQSAikuUUBCIiWU5BICKS5RQEIiJZLmuCYP2OWm6bs4KlG2tQ19siIu/pcg+UHanX1+/iF39fzc/+9jZD+xZywfEDmTVtJKUl+VGXJiISKetq/zquqKjwI32yeMe+Bp5etoU/L93CCyu3UVIQ52vnH8tVFUMwswxXKiLSeZjZq+5e0eqybAqCllZu3cPXHl7CgjU7OWVEX358zSQG9CrIQIUiIp3PoYIgtGsEZna3mVWZ2dI2lk83sxozWxi8bgmrltaM6d+D386ayq1XHM+yTbu5+o55bNy1vyNLEBHpFMK8WDwbmNHOOs+7+0nB69sh1tKqnBzj6inDuG/mKezY18DVv5zH+h21HV2GiEikQgsCd58L7Ajr9zNp0rA+3H/DB9hT18Q1d7ysMBCRrBL17aNTzWyRmf3ZzCa2tZKZzTKzSjOrrK6uDqWQE4b05v4bPsDe+ib+YfYCavY3hrIdEZHOJsogeA04xt1PBH4CPNLWiu5+h7tXuHtFeXmr4ypkxHGDe/GLj5/Mmm37+Oz9r9GYSIa2LRGRziKyIHD33e6+N/j8JJBrZmVR1dNs6qhSvnv58bywahu3PLpUD5+JSLcX2QNlZjYA2OrubmankAql7VHV09JHK4ayZts+fv7c24zp14NPnzEi6pJEREITWhCY2YPAdKDMzDYA/w7kArj7L4ArgX8ysyZgP3CNd6J/fn/xQ+N4a+tevvvkck4c2ouTj+kbdUkiIqHI2gfK0lGzv5GLfvICDU1Jnrj5DMrUHYWIdFGRPFDWHfQqzOXnH5vMjtoGPv+b10kku1ZoioikQ0HQjuMG9+I7l0zkxVXb+fEzK6MuR0Qk4xQEabh6yjAunzSYnzy7kgVrusQzciIiaVMQpOlbl0xkSJ8ivvCbhXrYTES6FQVBmnoU5PKja05iy+46bnm01X70RES6JAXBYZg8rA+fP2cMjy7cxMOvb4i6HBGRjFAQHKbPnj2aimP68O+PvsHW3XVRlyMictQUBIcplmPcdtWJNCSSfPWPS9QFhYh0eQqCIzCirJgvfXg8z66o4uHXN0ZdjojIUVEQHKHrTxtOxTF9+OZjaiISka5NQXCEYjnGf155AvVNSd1FJCJdmoLgKIwsL+Hmc8Yw542tvLhqW9TliIgcEQXBUZp5xgiG9i3k248vo0kD2YhIF6QgOEoFuTH+7YIJvLl1Dw++si7qckREDpuCIAM+PLE/p40q5YdPv8Wu2oaoyxEROSwKggwwM265aAK79zfyo7+qh1IR6VoUBBkyfkBPrp4ylAfmr2NLjW4nFZGuQ0GQQTdNH03CnTvmro66FBGRtCkIMmho3yIuPWkwD7yylm1766MuR0QkLQqCDLvp7FHUNyW564V3oi5FRCQtCoIMG1VewoXHD+TXL63RHUQi0iUoCELw2bNHs68hweyX1kRdiohIuxQEITh2YE/OPbY/97y4hn31TVGXIyJySAqCkNx09ihq9jfqaWMR6fRCCwIzu9vMqszskF1zmtkUM0uY2ZVh1RKFycP6cOrIvvzq+dXUNyWiLkdEpE1hnhHMBmYcagUziwG3AnNCrCMyN00fzdbd9TyiwWtEpBMLLQjcfS6wo53V/hn4A1AVVh1RmjamjImDevLLv68mkdSQliLSOUV2jcDMBgOXAb9IY91ZZlZpZpXV1dXhF5chZsZN00ezets+5ryxJepyRERaFeXF4h8B/+ru7Tagu/sd7l7h7hXl5eUdUFrmzDhuACPKirn9ubc10L2IdEpRBkEF8BszWwNcCfzczC6NsJ5QxHKMf5w2kiUba3h5dXstZSIiHS+yIHD3Ee4+3N2HA78HbnL3R6KqJ0yXTx5MaXEedz6vzuhEpPMJ8/bRB4F5wDgz22BmM83sRjO7MaxtdlYFuTE+fuoxPLOiilVVe6MuR0TkAPGwftjdrz2Mda8Pq47O4hNTj+H2v7/NXS+8w/cuPz7qckRE3qUniztIWUk+V0wezB9f28B2dVEtIp2IgqADzTxjBPVNSe57eW3UpYiIvEtB0IFG9+vB2ePKuW/eWuoa1e2EiHQOCoIONvOMkWzf18CTSzZHXYqICKAg6HCnjy5lZHkx985T85CIdA4Kgg5mZnxq6nAWrd/FwvW7oi5HRERBEIXLJw+mOC/Gr+etiboUEREFQRR6FORyxclDeGLRZt1KKiKRUxBE5JNTj6EhkeQ3C9ZHXYqIZDkFQURG9+vB6aNLuf/ltTQlklGXIyJZTEEQoU+cOpxNNXX8/a2uM8aCiHQ/CoIInXNsP8p75GuAexGJlIIgQrmxHK46eQjPrqhic83+qMsRkSylIIjYNVOGkXR4qHJD1KWISJZSEERsWGkRZ4wu47cL1muAexGJhIKgE7j2lGFs3LWf51fqorGIdDwFQSdw3oT+lBbn6aKxiERCQdAJ5MVzuPLkIfx1eRXVe/SksYh0LAVBJ3H55CEkks5f3tgSdSkikmUUBJ3E2P4ljCwv5snFGqdARDqWgqCTMDMuPH4g89/ZzjZ1RCciHUhB0IlccPxAkg5/WarmIRHpOAqCTmT8gB6MLCvmz0vVPCQiHafdIDCzq9KZ18o6d5tZlZktbWP5JWa22MwWmlmlmZ2RXsndl5lxwfEDmff2do1TICIdJp0zgq+mOe9gs4EZh1j+DHCiu58EfBq4M43f7Paam4fmvLE16lJEJEvE21pgZucDFwCDzezHLRb1BJra+2F3n2tmww+xfG+LyWJA/SsAxw7swYiyYp5cspnrPjAs6nJEJAsc6oxgE1AJ1AGvtng9Bnw4Exs3s8vMbAXwJ1JnBW2tNytoPqqsru7e3TCkmocGMG+1modEpGO0GQTuvsjd7wVGu/u9wefHgFXuvjMTG3f3h919PHAp8J1DrHeHu1e4e0V5eXkmNt2pnX/cQBJJ55kVVVGXIiJZIJ1rBE+bWU8z6wssAu4xs//KZBHuPhcYZWZlmfzdrmrioJ4M7l3IU7pOICIdIJ0g6OXuu4HLgXvc/WTg3KPdsJmNNjMLPk8G8oDtR/u73YGZcd6E/jy/sprahnYvx4iIHJV0giBuZgOBjwJPpPvDZvYgMA8YZ2YbzGymmd1oZjcGq1wBLDWzhcDPgKvdXReMAx+a0J/6piTPr9wWdSki0s21eddQC98G5gAvuvsCMxsJrGzvS+5+bTvLbwVuTavKLDRlRF96Feby1Btb+fDEAVGXIyLdWLtB4O4PAQ+1mF5N6l/zEqLcWA7njO/HMyu20pRIEo/pIXARCUc6TxYPMbOHg6eEt5rZH8xsSEcUl+3Om9CfXbWNLFiTkZu0RERalc4/M+8hddvoIGAw8HgwT0J25thy8uI5PLVMndCJSHjSCYJyd7/H3ZuC12yg+9/M3wkU58eZNrqMp5dtRdfRRSQs6QTBNjP7uJnFgtfH0W2eHea8Cf3ZsHM/yzbvjroUEemm0gmCT5O6dXQLsBm4kkN0ByGZdd6E/sRyjCeXqGtqEQlHu0Hg7uvc/WJ3L3f3fu5+qbuv7YjiBEpL8jltVCmPL9qs5iERCUU6dw3da2a9W0z3MbO7wy1LWrroxEGs21HL4g01UZciIt1QOk1DJ7j7ruaJoMO5SeGVJAf78MQB5MaMxxdtiroUEemG0gmCHDPr0zwRdD6XzhPJkiG9CnM5a2w/nli8mWRSzUMiklnpBMEPgZfM7Dtm9m3gJeA/wy1LDnbRiQPZsruOyrV6uExEMiudi8W/JtWlxFagGrjc3e8LuzA50LnH9qcgN0fNQyKScWl1YOPuy9z9p+7+E3dfFnZR8n7F+XHOObY/Ty7ZTFMiGXU5ItKNqCezLuSiEwaxfV8DL72t5/lEJHMUBF3I9HHllOTH+dNiPVwmIpmTVhCY2TFmdm7wudDMeoRblrSmIDfGucf2Y86yLTSqeUhEMiSdB8r+Efg98Mtg1hDgkTCLkrZdcPxAdtU2qnlIRDImnTOCzwKnA7sB3H0l0C/MoqRtZ45NNQ89qeYhEcmQdIKg3t0bmifMLA7oqaaIqHlIRDItnSD4u5l9DSg0s/NIDVv5eLhlyaFceMIgNQ+JSMakEwRfIfUg2RLgM8CTwNfDLEoObdqYMjUPiUjGpDN4fRL4VfCSTqAgN8Z5E/ozZ9kW/l/iOHI1sL2IHIV07hpaYmaLD3o9b2b/bWalHVGkvF/z3UMvrtoWdSki0sWl04von4EE8EAwfU3wvhuYDVyU+bKkPWeOLaNHQZzHF21m+jjdxCUiRy6dNoXT3f2r7r4keP0bMN3dbwWGt/UlM7vbzKrMbGkbyz/W4gzjJTM78ch2ITvlx2PMmDiAOW9soa4xEXU5ItKFpRMEJWb2geYJMzsFKAkmmw7xvdnAjEMsfwc4y91PAL4D3JFGLdLCJScNZm99E39bURV1KSLShaXTNHQDcLeZlQBGqknoBjMrBr7X1pfcfa6ZDT/E8pdaTL5M6ollOQxTR5VSVpLPY4s2cf7xA6MuR0S6qHTuGloAHG9mvQBrOWwl8LsM1TGT1LWIVpnZLGAWwLBhwzK0ya4vlmN85ISBPPDKOnbXNdKzIDfqkkSkC0q307kLST1DcLOZ3WJmt2SqADM7m1QQ/Gtb67j7He5e4e4V5eXlmdp0t3DxSYNoaEoyZ+mWqEsRkS4qndtHfwFcDfwzqaahq4BjMrFxMzsBuBO4xN31mOwRmDS0N0P7FvKYRi4TkSOUzhnBae7+SWCnu38LmAoMPdoNm9kw4I/AJ9z9raP9vWxlZlx84iBeXLWN6j31UZcjIl1QOkFQF7zXmtkgoBEY0d6XzOxBYB4wzsw2mNlMM7vRzG4MVrkFKAV+bmYLzazyCOoXUncPJR3+tFhnBSJy+NK5a+hxM+sN3Aa8Rqrn0Xa7m3D3a9tZfgOpO5LkKI3t34NjB/bk4YWbuP70djNaROQAhzwjMLMc4Bl33+XufyB1bWC8u2fsYrFkxmWTBrFo/S5WV++NuhQR6WIOGQRBh3M/bDFd7+41oVclh+2SkwZjBo8sVPOQiByedK4RPGVmV5iZhV6NHLH+PQs4fVQZj7y+EXeNGyQi6UsnCP6F1GA0DWa228z2mNnukOuSI3DppMGs21HLa+t2Rl2KiHQh7QaBu/dw9xx3z3X3nsF0z44oTg7PjOMGUJCbwx9f2xh1KSLShaTzQJmZ2cfN7BvB9NCg4znpZEry43xowgCeWLyZhiaNZywi6UmnaejnpB4iuy6Y3gv8LLSK5KhcNnkwNfsbee5N9UgqIulJJwg+4O6fJXiwzN13AnmhViVHbNroMspK8vjDaxuiLkVEuoh0gqDRzGKkHiTDzMoBtTt0UvFYDpdNGswzy6vYvlddTohI+9IJgh8DDwP9zOw/gBeA74ZalRyVqyqG0pR0PVMgImlJ566h+4EvkxqEZjNwqbs/FHZhcuTG9u/BiUN68VDlej1TICLtSueuof8B+rr7z9z9p+6+vAPqkqN0ZcVQVmzZwxub9MiHiBxaOk1DrwFfN7NVZnabmVWEXZQcvYtPGERePIeHKtdHXYqIdHLpNA3d6+4XAKcAbwG3mtnK0CuTo9KrKJcPTxzAIws3UdeYiLocEenE0hqqMjAaGA8MB1aEUo1k1FUnD6FmfyN/Xb416lJEpBNL5xpB8xnAt4E3gJPd/aLQK5OjdvroMgb1KuB3lXqmQETals7ANO8AU919W9jFSGbFcowrK4byk2dXsnHXfgb3Loy6JBHphNK5RvALIGFmp5jZmc2vDqhNMuCjFUMA+N0CXTQWkdal0zR0AzAXmAN8K3j/ZrhlSaYM6VPEtDHlPFS5nkRSzxSIyPulc7H488AUYK27nw1MAqpDrUoy6popQ9lUU8fzK3XYROT90gmCOnevAzCzfHdfAYwLtyzJpHOP7U/f4jx+q+YhEWlFOkGwwcx6A48AT5vZo4A6selC8uI5XDF5ME8v20r1HnVEJyIHSudi8WXuvsvdvwl8A7gLuDTswiSzrp6S6ojuj+qeWkQOcjgPlOHuf3f3x9y9ob11zexuM6sys6VtLB9vZvPMrN7Mvng4dcjhG92vB1OG9+F/56/VRWMROcBhBcFhmg3MOMTyHcDNwA9CrEFa+PTpI1i/Yz9PL9sSdSki0omEFgTuPpfUH/u2lle5+wKgMawa5EAfmjiAoX0L+dXz70Rdioh0ImGeEWSMmc0ys0ozq6yu1i2QRyqWY3z69BG8unYnr63bGXU5ItJJdIkgcPc73L3C3SvKy8ujLqdLu6piKD0K4tylswIRCXSJIJDMKcmPc90HhvHnpZtZv6M26nJEpBNQEGSh608bTo4Z97y4JupSRKQTCC0IzOxBYB4wzsw2mNlMM7vRzG4Mlg8wsw3Av5AaAW2DmfUMqx55z8BehVx84iAeeGUtVbvroi5HRCKWTjfUR8Tdr21n+RZgSFjbl0O7+ZwxPLZoEz95dhXfufS4qMsRkQipaShLDS8r5qNThvLgK+tYt13XCkSymYIgi938wTHEcowf/fWtqEsRkQgpCLLYgF4FXH/acB5euJE3t+yJuhwRiYiCIMvdeNYoSvLi3DbnzahLEZGIKAiyXJ/iPG6cPoq/Lt/K396sirocEYmAgkC4YdoIRpUXc8ujS6lrTERdjoh0MAWBkB+P8Z1Lj2P9jv389NlVUZcjIh1MQSAAnDaqjMsnDeaXc99mVdXeqMsRkQ6kIJB3fe3CYynMjfH1R5bgrsFrRLKFgkDeVVaSz1fOP5aXV+/g/vnroi5HRDqIgkAOcO0pQ5k2pozvPrlcTxyLZAkFgRzAzLj1ihOImfGl3y8iqfGNRbo9BYG8z6DehXzjIxOY/84Ofj1vTdTliEjIFATSqqsqhnD2uHK+/5cVrKpS9xMi3ZmCQFrV3ERUlBfncw+8rgfNRLoxBYG0qV/PAn541Yms2LKH7z25POpyRCQkCgI5pLPH92PmGSO4d95anl62NepyRCQECgJp15dnjGPioJ586feLNOC9SDekIJB25cdj/PS6ySSSzqz7XqW2oSnqkkQkgxQEkpYRZcX85NpJvLllN196aLG6oBDpRhQEkrbp4/rxrzPG86clm/n5c29HXY6IZEg86gKka5l15kiWbd7ND556k1HlJcw4bkDUJYnIUdIZgRyW5ucLThzSmy/89nUWrt8VdUkicpQUBHLYCnJj3PmpCsp75HPDvQt0J5FIFxdaEJjZ3WZWZWZL21huZvZjM1tlZovNbHJYtUjmlZXkc8/1p9CYcP5h9gJ27muIuiQROUJhnhHMBmYcYvn5wJjgNQu4PcRaJASj+5Xwy0+czLodtXzszvkKA5EuKrQgcPe5wI5DrHIJ8GtPeRnobWYDw6pHwnHqyFJ+9ckK3q7ey7W/epnte+ujLklEDlOU1wgGA+tbTG8I5r2Pmc0ys0ozq6yuru6Q4iR9Z40t565PTeGdbfu47lfzqdpTF3VJInIYogwCa2Veq08pufsd7l7h7hXl5eUhlyVH4owxZdxz/RTW7ajl0p++yPLNu6MuSUTSFGUQbACGtpgeAmyKqBbJgNNGl/HQjVNJOlx5+0s8u0Kd1Il0BVEGwWPAJ4O7h04Fatx9c4T1SAYcN7gXj37udEaWlzDz3kp+9rdVGu5SpJML8/bRB4F5wDgz22BmM83sRjO7MVjlSWA1sAr4FXBTWLVIx+rfs4DffWYqF50wiNvmvMn1sxfoIrJIJ2ZdrfOwiooKr6ysjLoMSYO78+Ar6/nm42/QpyiX//7oSZw2uizqskSykpm96u4VrS3Tk8USGjPjug8M4+GbTqM4L851d87n648sYV+9urEW6UwUBBK6iYN68aebpzHzjBHcP38dH/7RXJ5fqduARToLBYF0iMK8GN/4yAQe+sxUcmM5fOKuV/jCb15nm64diEROQSAdqmJ4X/78+WncfM4Y/rRkMx/8wXPc9/JamhLJqEsTyVoKAulwBbkx/uW8sfz582cyYVBPvvHIUi788Qu8sHJb1KWJZCUFgURmdL8SHvzHU7n9Y5OpbWzi43fN5/p7XmHJhpqoSxPJKrp9VDqFusYEs19aw+3PvU3N/kbOm9CfL5w7homDekVdmki3cKjbRxUE0qnsqWvk7hfWcOcLq9lT18Tpo0u5YdpIpo8tx6y17qlEJB0KAulyavY38sD8dcx+6R227q5nVHkx154yjCsmD6FPcV7U5Yl0OQoC6bIampI8sXgT//vyWl5bt4u8WA4fmtifyyYNZtqYcvLiuswlkg4FgXQLK7bs5jevrOfRhRvZWdtI76Jczj9uAB+aOIDTRpWSH49FXaJIp6UgkG6lMZHkhZXbeHThRp5etpV9DQlK8uOcNbacs8aWc+bYcgb0Koi6TJFO5VBBEO/oYkSOVm4sh7PH9+Ps8f2ob0rw0tvbeeqNrTyzfCt/WpLqyXxs/xJOHVnKqSNLOWVEX8pK8iOuWqTz0hmBdBvuzptb9zD3rWqeX7mNV9fupLYhAcDw0iImH9OHycP6cMKQXowb0ENNSZJV1DQkWakxkWTJxhrmr97Ba+t28vq6nWzb2wBAbswYN6AHEwb2ZPyAnhw7sCdj+5dQqjMH6abUNCRZKTeWw+RhqbMASJ0xrN+xnyUba1iysYalG2t4ZnkVv6vc8O53+hbnMbq8hFH9ihlRVsyIshKOKS1iWN8iCnJ1BiHdk4JAsoaZMay0iGGlRVx4wsB351fvqWf55t2srNrLqqo9rNy6l6fe2Mr2fQ0HfL9/z3yG9ilicJ9ChvQpZFDvQgb1KmRg7wIG9iqkZ0FcD71Jl6QgkKxX3iOf8h6pu41aqqltZPW2vazbUcu67bWs3VHLhp21vLp2J08s3kzioLGYC3JzGNCzgH49CijvmU95SX7qt0vyKS3Jo7Qkn9LiPPoW51GUF1NoSKehIBBpQ6+iXCYN68OkoGmppaZEkqo99Wyu2c+mXXVsqalj6+46tuyuo2p3Pcs27aZ6Tz172xiNLT+eQ5+iPPoU59GnKJc+RXn0LMyld1EuvQpz6VkQvBfG6VmQS4+COD2C9/x4jkJEMkpBIHIE4rGcVNNQ70JOPqbt9fY3JNi2t55te+vZvreBHbUN7NiXeu3c18DO2gZ21jayYstuavY3sau2gabkoW/gyI0ZJflxSgrilOTnUpwXozg/Tkl+nKLgc3F+jKK8OIW5MYrzYxTmxSnKjVGYF7xyYxQF7/m5qffcmClgspSCQCREhXkxhvYtYmjforTWd3f2Nyao2d9Izf5G9tQ1saeukd37U+976pvYU9fE3rom9tU3sac+9b6ztoH1O2vZ35BgX30T+xoS72u6ak+OQWFujILglZ+bQ348RkFuDvnxnNS8Fu/58Rh58dSyvOZXLIf83Bj5sQPnNX/OjaXWz43lkBuzA6fjwbycHHJyFEgdSUEg0omYGUV5cYry4gzsVXjEv+PuNCSS1NYn2NfQRF1jgtqG1Gt/Q4K6xgT7g3l1je9N1zUmg+kk9U3vvdc3Jtmxr4G6xgQNTUnqg1fqc4LGRGZvQ4/npEIiHjPygvdUeKTCIp7zXpA0L4vlvDc/HsshN8eIt/gcC5bFclrMC4InFqwbyzlwOt78OVg3nmPvbufd+c3ftfeW5eRwwHvM3vt+jgXvnSjsFAQi3ZCZkR+PkR+PdUhvrYmk05hoDohUWDQ0JWlIpN6blzUm/N1lTcn31mlsXpZ4b1nzuk3JJI1Nqd9vTDpNidTvNSb83WV7m5qCGlLLm4J6mprXSfi7NTYl/bDPlsISC4Ik1iIcmt+bg6X5lWNw7SnDuGHayIzXEWoQmNkM4H+AGHCnu3//oOXHAHcD5cAO4OPuvuF9PyQinVrqj1UseNYiN+py2uWeCoOm4JVIOI3J5HvzEqnPzeGSmp98dzrpB66X9Bbzg/UT3vzbSRIOyWRqG8mkk0hCIpkk4S0+B+9Nwe8lgvWat5VMemhdpYQWBGYWA34GnAdsABaY2WPuvqzFaj8Afu3u95rZB4HvAZ8IqyYREUidMaWajaKupHMIszP3U4BV7r7a3RuA3wCXHLTOBOCZ4PPfWlkuIiIhCzMIBgPrW0xvCOa1tAi4Ivh8GdDDzEoP/iEzm2VmlWZWWV1dHUqxIiLZKswgaO2S+MFXaL4InGVmrwNnARuB9z2B4+53uHuFu1eUl5cfvFhERI5CmBeLNwBDW0wPATa1XMHdNwGXA5hZCXCFu9eEWJOIiBwkzDOCBcAYMxthZnnANcBjLVcwszIza67hq6TuIBIRkQ4UWhC4exPwOWAOsBz4nbu/YWbfNrOLg9WmA2+a2VtAf+A/wqpHRERap4FpRESywKEGpgmzaUhERLqALndGYGbVwNoj/HoZsC2D5XQV2bjf2bjPkJ37nY37DIe/38e4e6u3XXa5IDgaZlbZ1qlRd5aN+52N+wzZud/ZuM+Q2f1W05CISJZTEIiIZLlsC4I7oi4gItm439m4z5Cd+52N+wwZ3O+sukYgIiLvl21nBCIichAFgYhIlsuaIDCzGWb2ppmtMrOvRF1PGMxsqJn9zcyWm9kbZvb5YH5fM3vazFYG732irjUMZhYzs9fN7IlgeoSZzQ/2+7dBn1fdhpn1NrPfm9mK4JhPzYZjbWb/J/jve6mZPWhmBd3xWJvZ3WZWZWZLW8xr9fhayo+Dv2+LzWzy4WwrK55PWKwAAATLSURBVIKgxWhp55MaDOdaM5sQbVWhaAL+r7sfC5wKfDbYz68Az7j7GFIDAXXLIAQ+T6pfq2a3Av8d7PdOYGYkVYXnf4C/uPt44ERS+96tj7WZDQZuBirc/ThSw+BeQ/c81rOBGQfNa+v4ng+MCV6zgNsPZ0NZEQSkN1pal+fum939teDzHlJ/GAaT2td7g9XuBS6NpsLwmNkQ4ELgzmDagA8Cvw9W6Vb7bWY9gTOBuwDcvcHdd5EFx5pU9/mFZhYHioDNdMNj7e5zSY3l3lJbx/cSUsP+uru/DPQ2s4HpbitbgiCd0dK6FTMbDkwC5gP93X0zpMIC6BddZaH5EfBlIBlMlwK7gl5wofsd85FANXBP0Bx2p5kV082PtbtvJDXW+TpSAVADvEr3PtYttXV8j+pvXLYEQTqjpXUbwSA/fwC+4O67o64nbGb2EaDK3V9tObuVVbvTMY8Dk4Hb3X0SsI9u1gzUmqBN/BJgBDAIKCbVLHKw7nSs03FU/71nSxC0O1pad2FmuaRC4H53/2Mwe2vzaWLwXhVVfSE5HbjYzNaQavb7IKkzhN5B8wF0v2O+Adjg7vOD6d+TCobufqzPBd5x92p3bwT+CJxG9z7WLbV1fI/qb1y2BEG7o6V1B0G7+F3Acnf/rxaLHgM+FXz+FPBoR9cWJnf/qrsPcffhpI7ts+7+MeBvwJXBat1qv919C7DezMYFs84BltHNjzWpJqFTzawo+O+9eb+77bE+SFvH9zHgk8HdQ6cCNc1NSGlx96x4ARcAbwFvA/8WdT0h7eMZpE4HFwMLg9cFpNrLnwFWBu99o641xP8NpgNPBJ9HAq8Aq4CHgPyo68vwvp4EVAbH+xGgTzYca+BbwApgKXAfkN8djzXwIKnrII2k/sU/s63jS6pp6GfB37clpO6qSntb6mJCRCTLZUvTkIiItEFBICKS5RQEIiJZTkEgIpLlFAQiIllOQSBZy8xeCt6Hm9l1Gf7tr7W2LZHOSLePStYzs+nAF939I4fxnZi7Jw6xfK+7l2SiPpGw6YxAspaZ7Q0+fh+YZmYLg77uY2Z2m5ktCPp2/0yw/vRgvIcHSD20g5k9YmavBv3jzwrmfZ9U75gLzez+ltsKnvy8LehLf4mZXd3it59rMb7A/cGTsyKhi7e/iki39xVanBEEf9Br3H2KmeUDL5rZU8G6pwDHufs7wfSn3X2HmRUCC8zsD+7+FTP7nLuf1Mq2Lif1RPCJQFnwnbnBsknARFJ9xLxIqg+lFzK/uyIH0hmByPt9iFS/LQtJdeNdSmrAD4BXWoQAwM1mtgh4mVSnX2M4tDOAB9094e5bgb8DU1r89gZ3T5LqHmR4RvZGpB06IxB5PwP+2d3nHDAzdS1h30HT5wJT3b3WzJ4DCtL47bbUt/icQP//lA6iMwIR2AP0aDE9B/inoEtvzGxsMOjLwXoBO4MQGE9qeNBmjc3fP8hc4OrgOkQ5qVHGXsnIXogcIf2LQyTVe2dT0MQzm9RYwMOB14ILttW0PvThX4AbzWwx8Cap5qFmdwCLzew1T3WJ3exhYCqwiFRPsV929y1BkIhEQrePiohkOTUNiYhkOQWBiEiWUxCIiGQ5BYGISJZTEIiIZDkFgYhIllMQiIhkuf8PG7THG+H1kUUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(results)\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('average cost')\n",
    "plt.savefig(\"adp.pdf\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
